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Abstract 

Bayesian estimation strategies represent the most fundamental formulation of the state estimation problem available, and apply 
readily to nonlinear systems with non-Gaussian uncertainties. The present paper introduces a novel method for implementing 
grid-based Bayesian estimation which largely sidesteps the severe computational expense that has prevented the widespread 
use of such methods. The method represents the evolution of the probability density function (PDF) in phase space, p x (x',t), 
discretized on a fixed Cartesian grid over all of phase space, and consists of two main steps: (i) Between measurement times, 
p x (x',i) is evolved via numerical discretization of the Kolmogorov forward equation, using a Godunov method with second- 
order corner transport upwind correction and a total variation diminishing flux limiter; (ii) at measurement times, p x (x', t) is 
updated via Bayes' theorem. Computational economy is achieved by exploiting the localised nature of p x (x', t). An ordered list 
of cells with non-negligible probability, as well as their immediate neighbours, is created and updated, and the PDF evolution 
is tracked only on these active cells. 
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1 Introduction 

Bayesian estimation strategies are the most general class 
of solutions to the state estimation problem, and apply 
readily to nonlinear systems where information about 
the state is represented by a probability density func- 
tion (PDF) of general form. In this paper we introduce 
a novel, computationally cheap method for implement- 
ing grid-based Bayesian estimation that exploits the fact 
that the PDF is usually negligible in most of phase space, 
while avoiding many of the disadvantages of other meth- 
ods. The idea of grid-based Bayesian estimation dates 
back at least to Stratonovich (1959,1960). The equa- 
tions underlying the algorithm are laid out clearly in 
Jazwinski (1970, p. 164), and are summarised below. 
However, numerical implementation of these equations 
has only been attempted sporadically in the half cen- 
tury since, for instance by Kramer et al. (1988), Terwi- 
esch & Agarwal (1994) and Ungarala et al. (2006). Grid- 
based Bayesian methods typically suffer from the twin 
burdens of high computational cost and a finite domain 
size; indeed, Arulampalam et al. (2002), in their other- 
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wise insightful review of particle filter methods, all but 
dismiss grid-based methods in §IV.B. We believe that 
this level of pessimism on this class of methods is un- 
warranted. The algorithm developed, dubbed GBEES 
(Grid-based Bayesian Estimation Exploiting Sparsity), 
provides a means of efficient computation by building on 
an accurate integration scheme for hyperbolic systems, 
and a novel gridding scheme over all of phase space. 

2 Grid-based Bayesian estimation exploiting 
sparsity 

Consider the state estimation of the nonlinear system 
-=f(x,w), y = h(x,v). (1) 

The grid-based Bayesian estimation method is best vi- 
sualized as an evolution of the PDF of the state estimate 
x discretized on a fixed grid over all of phase space R"; 
assuming the state x develops according to the nonlin- 
ear equation (1), the method consists of two relatively 
straightforward steps (for details, see Jazwinski 1970, 
p. 164): 

(i) Between measurement times, the PDF itself, p x (x', t) , 
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is marched via discretization of the Kolmogorov forward 
equation (also called the Fokker-Planck equation) 

dp x (x',f) _ dfj(x',t)p x (xL',t) 1 d 2 qijp x (x',t) 
dt ~ dx\ 2 dx'tdx'j ' 

where summation over repeated indices is implied and 
Qij is the (i, j)th element of the spectral density, Q, of the 
state disturbances (note that, in the special case that Q 
is diagonal and the state disturbance is independent of 
x, this is just a diffusion term). Risken (2002) reviews a 
number of methods for solving this equation, including 
analytic methods for special cases and eigenfunction ex- 
pansions (focusing on the stationary solution) , however 
state-space and time discretisation for a non-stationary 
solution (as performed here) is only mentioned briefly, 
and assumed to apply only to finite domains. An ac- 
curate numerical method for marching this equation in 
time in the case that Q = is outlined in §3; adding an 
appropriate term to this discretization to apply diffusion 
to the PDF (to account for Gaussian state disturbances) 
is straightforward, as discussed in §3.2. 

(ii) At the measurement times tk, the PDF is updated 
via Bayes' theorem (Bayes, 1763), 

„ fv > f n Py(vfc|x')Px(x , ,t fc -) , . 

p x (x,t k +) = ^ , (3) 

where p x (x',t k +) denotes the a posteriori PDF (after 
accounting for the measurement y^), p y (vfc|x') denotes 
the uncertainty associated with the measurement (which 
may or may not be Gaussian in x'), p x (x', tk- ) denotes 
the a priori PDF (before accounting for the measure- 
ment yfc), and C is an appropriate normalization con- 
stant, which is selected for every measurement update 
to normalize the discretization of p x (x',tj.+ ) such that 
its integral over phase space is unity. 

To understand how the continuous state-space is dis- 
cretized, recall that the cumulative distribution function 
(CDF) of a random real vector x, denoted / x (x), maps 
x e R™ to the real interval [0, 1] that monotonically in- 
creases in each of the components of x, and is defined 

/ x (x) = P(x 1 <x 1 ,x 2 <x 2 ,...,x n < x n ), 

where x is some particular value of the random vector 
x and P(S) denotes a probability measure that the con- 
ditions stated in S are true. For any random vector x 
whose CDF is differcntiable everywhere, the probability 
density function (PDF) p x (x') > is a scalar function 
of x' defined such that 

rui rn 2 r* n 
/x(x) = / / •••/ p^)dx' 1 dx' 2 ■■■dx' n , 



dxi dx 2 ' ' ' 9x,, 



For small |Ax'|, the quantity p x (x')Aa;' 1 Ax' 2 • ■ • Ax' n 
represents the probability that the random vector x 
takes some value within a small rectangular region cen- 
tered at the particular value x' and of width Ax- in 
each coordinate direction e, . 

The method we have developed maintains a list of active 
cells on the grid over all of phase space in order to limit 
both the computational effort and the memory storage 
required in the numerical simulation. This list includes 
all cells in the discretization for which the PDF is greater 
than a given threshold, as well as all cells which, though 
they may or may not themselves exceed this threshold, 
are either one of the two immediate neighbor cells, in 
each of the n coordinate directions, of those cells which 
exceed the threshold, or are one of the four neighbor 
cells, in each of the n C 2 pairs of coordinate directions, 
which touch a corner of the cells which exceed the thresh- 
old. At each time step, cells are added to and removed 
from this list as appropriate, and the fluxes initialized 
and updated on every interior boundary between adja- 
cent cells in the list. When performing a computation 
restricted to an evolving list of active grid cells of this 
sort, the relative position of the various cells in the list 
is needed frequently. This may be determined efficiently 
by keeping in each list record a pointer to the two imme- 
diate neighbor cells in each coordinate direction in the 
list, if these neighbor cells are present in the list, or to 
a null record if not, and updating these pointers appro- 
priately as records are added to and removed from the 
list 1 . These pointers interconnecting the list facilitate 
rapid computation of the numerical discretization given 
in (5) in the next section. The most expensive step in 
maintaining this list of neighbor cells is making the ap- 
propriate connections when a new record is added to the 
list. Though this may be accomplished by scanning the 
entire list, this approach becomes prohibitively expen- 
sive as the length of the list grows to thousands of cells. 
Instead, we keep the list ordered by its indices (e.g., in 
a phase space with n = 3, ordered first by i, then by j, 
then by k) , and store the elements of the list as a binary 
tree. This allows the time-limiting search step to proceed 
at O(ATlogiV) operations, where N is the number of list 
elements. Conveniently, this list ordering and searching 
can be handled using the C++ Standard Template Li- 
brary map container. 

3 Accurate numerical integration of the Kol- 
mogorov forward equation 

The PDE governing the evolution of the PDF in the 
present problem is given by (2). If Q = 0, the equation 
is hyperbolic; if Q > 0, the equation, strictly speak- 
ing, changes type to elliptic. In practice, however, Q is 
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usually relatively small. It is thus fitting to design a 
numerical method for accurate simulation of (2) based 
on a proven algorithm for accurate simulation of hyper- 
bolic PDEs. Fortunately, the fluid mechanics commu- 
nity has focused on the development of high performance 
computing techniques for numerical simulation of such 
"convection-dominated" problems for over 40 years, and 
these techniques are now quite refined and well under- 
stood. The numerical method best suited to the present 
problem is somewhat involved; a comprehensive review 
of this class of methods is given in LeVeque (2002). To 
focus this discussion, consider first the two-dimensional, 
linear, hyperbolic PDE in conservation form 

dpjx, y, t) ^ du(x,y)p(x,y,t) _ dv(x,y)p(x,y,t) 
dt dx dy 

(4) 

noticing that higher-dimensional cases follow as an ob- 
vious extension. Following Chapters 4, 6, 9, 19, and 20 
of LeVeque (2002), we implement a Godunov-type finite 
volume method by writing (4) on a uniform Cartesian 
2D mesh (with constant Ax and Ay) in the form 
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1/2, j are determined, 

for all i and j, by first initializing 
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where u + = max(u, 0), u~ = min(u, 0), etc, then ap- 
plying the corner transport upwind (CTU) terms by up- 
dating, for all i and j, 
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and the flux limiter function <p(9) e [0, 2] is selected 
as one of several possible choices, including the raono- 
tonized central- difference (MC) limiter and the van Leer 
limiter: 



MO. <t>{6) = max{0, min[(l + 9) /2, 2, 29)}, 

van Leer: <j)(9) = (9 + \9\) / (1 + \9\) . 

Note that exact conservation of the discrete approxima- 
tion of the integral of p over phase space, as implied by 
the continuous formulation in (4), follows immediately 
from (5). 
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3. 1 Numerical analysis 

In regions characterized by smooth variation of p, 9 rts 
1 and (f)(9) rts 1, and the algorithm described in §3 
is amenable to straightforward numerical analysis. For 
simplicity, consider here the ID test problem 
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di 



dp 
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where u is a positive or negative constant. In this case, 
the discretization described above reduces to 
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and thus 
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Now applying to this equation the multidimensional 
Taylor series expansion, 
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and rearranging appropriately, gives 



(p i )r = -^x)r-^*(p«x n +^ 

+ C((At) 2 ,(Ax) 2 ,Aa;Ai). 



Differentiating (6) with respect to t and inserting (6) 
into the RHS of the result, it is seen that the second 
and third terms on the RHS of the above expression 
cancel. Thus, in regions of smooth variation of p, the 
proposed scheme is second- order accurate in both space 
and time. 2 A similar analysis follows for problems in 
higher dimensions. 



3.2 Accounting for diffusion 



A diffusion term is easily added to the discretization 
given in (5) in a second-order central finite difference 
fashion simply by updating the fluxes such that, for all 
i and j , 



Aw™_ 1/9 . 
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where /i is the coefficient of the diffusion term that is ap- 
plied numerically. The flux limitcr functions mentioned 
at the end of §3 are designed to reduce the algorithm, 
locally, to a first-order spatial behavior while applying 
sufficient numerical diffusion in regions of large local 
curvature of p on the grid, to provide a total variation 
diminishing (TVD) solution (that is, preventing spuri- 
ous oscillations with new local minima and maxima). 
We may compensate for the diffusion introduced by the 
numerical discretization of the convective terms simply 
by appropriately reducing the diffusion /2 applied in the 
numerical simulation of (2). 
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Fig. 1. Computation of a sensitive test system (governed by 
(4) with u(x, y) = y and v(x, y) — —x), using the method in 
§3, taking Q = 0, Ax = 0.01, At = 0.001, and a threshold 
of 10 (that is, tracking numerically only those cells with 
p > 0.001 and their immediate neighbors). The distributions 
are compared in cross section in the plane x = 0, with the 
black solid line denoting the exact solution at t — 2tt for 
Q — 0, the blue dashed line denoting the numerical solution 
at t = 2n using Q = 0, Ax = 0.01, At = 0.001 and the 
red dot-dashed line denoting the exact solution at t — 2ir 
for Q = 8 x 10~ 5 . The leading-order error in the numerical 
solution at t — 2n is a small amount of diffusion of p. 

3.3 Validation 

A simple yet sensitive numerical test of the algorithm 
is given in Figure 1; this numerical test was taken with 
u = y and v = —x in order to give simple solid body rota- 
tion about the origin, as suggested by LeVeque (2002). If 
Q = 0, the exact solution of the test problem considered 
in Figure 1 at t = 2tt, after a single rotation of the sys- 
tem about the origin, is simply the initial condition. For 
the case Q = 2/iI where \x is constant and positive, the 
exact solution of this problem at t = 2-7T may be obtained 
analytically by means of Fourier transforms. As seen by 
comparing Figure 1 to Figure 20.5 of LeVeque (2002), 
the result obtained via the GBEES approach is essen- 
tially identical to that obtained using the complete grid 
when sufficiently small threshold, time step and state- 
space discretization is used. The information loss due 
to the discretization scheme may be quantified via the 
Kullback-Liebler divergence (Kullback & Liebler, 1951), 
Dkl{Pqi Pe)t where a distribution P e is used to approx- 
imate the true distribution P . The Kullback-Liebler di- 
vergence using the simulation in Figure 1 when com- 
pared to the true analytic solution is 0.089 bits (with the 
distributions normalised to integrate to unity) , whereas 
the divergence from the true case to the diffusion case 
calculated analytically for Q — 8 x 10~ 5 is 0.085 bits 3 . 

2 Meaning that the error is bounded by a term proportional 
to {Ax) 2 in space and (At) 2 in time, giving convergence of 
0(1/N 2 ). 

3 This value of Q was chosen so that the divergence from 
the true case to the diffusion case was close to that of the 
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Fig. 2. The rapid transformation of a PDF from Gaussian 
to highly non-Gaussian in the Lorenz system, with no mea- 
surement updates. Visualized are p = 0.005, p = 0.0005, and 
p = 0.00005 isosurfaces of the PDF in phase space at t = 
(in the upper-right), t = 0.2, t = 0.4, t = 0.6, t = 0.8, and 
t = 1. The simulation was performed with a variable time 
step At < 0.001, a grid spacing of Ax — 0.25, and a threshold 
of e = 10 -6 . This simulation required less than 40 seconds of 
computation on a 2009 vintage Apple laptop computer (2.4 
GHz Intel Core 2 Duo) using a single-threaded CH — h imple- 
mentation of the present algorithm tracking about 50,000 
active cells (of 5.3 x 10 6 total cells in the domain shown) at 
t = 1. The modest memory requirements are proportional 
to the number of cells. 

The divergence from the true solution to a simulation 
using a truncation threshold of 10 -16 (not shown) is al- 
most the same as the more aggressively truncated ex- 
ample and is visually indistinguishable. As evident by 
comparing the numerical solution at t = 2?t in Figure 1 
to the initial condition, the discretization described in 
§3 introduces a small numerical error in regions of high 
curvature. However, by comparing the numerical solu- 
tion, for /J, = 0, Ay = 0.001 and Ax = 0.01, to the exact 
solution, for // = 4 x 10 -5 , it is evident that the leading- 
order error of the numerical discretization is just a bit 
of additional diffusion, the level of which may be deter- 
mined by a suitable minimisation process. 

4 Numerical results 

A Bayesian approach is justified when the uncertainty of 
the estimate is significantly non-Gaussian, such as in the 
estimation of a nonlinear system with relatively large un- 
certainty, leaving us with particle filtering or grid-based 
methods; what is perhaps still uncertain is the numerical 



true case to the numerical solution. 



Fig. 3. The same simulation as Figure 2 but with measure- 
ments of X3 (vertical axis), with Gaussian uncertainty, at 
every time step. The black line is the 'true' state which gen- 
erates the measurements. Note that, by t = 1, the PDF 
splits into two concentrated regions. This simulation required 
about 4 seconds of computation with the same hardware, 
tracking about 4,000 active cells (of 5.3 x 10 6 total cells in 
the domain shown) at t — 1. 

tractability of a grid-based approach when one exploits 
the sparsity of the PDF in the manner described in §2. 
Thus, in order to test the efficiency of the GBEES algo- 
rithm, as well as to demonstrate how it can capture with 
unprecedented accuracy the evolution of a non-Gaussian 
PDF, we have applied the GBEES algorithm to the es- 
timation of the three-state Lorenz system 
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with a = 4, b = 1, and r — 48. For these parameter val- 
ues, the system is chaotic, and the attractor takes the 
familiar form indicated by the green line in Figure 2. 
Also illustrated in this Figure 2 is the evolution of an 
initially (at t = 0) Gaussian PDF p x (x',£), the evolu- 
tion of which is governed by the Kolmogorov equation 
(2), with no measurement updates applied and no added 
process noise (diffusion). The distribution narrows sig- 
nificantly in the direction normal to the attractor, and 
spreads out rapidly in the direction of the maximum lo- 
cal Lyapunov exponent along the attractor; by t = 1, 
the PDF is highly non-Gaussian. Note also in the t = 0.8 
and t = 1 isosurfaces the remarkable division of the PDF 
into two distinct lobes in the vicinity of the x 3 axis (the 
vertical coordinate axis in the figures), which is invari- 
ant and unstable in the Lorenz system. 
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Figure 3 represents the evolution of the PDF when 
measurements (with Gaussian uncertainty) of X3 arc 
taken at every time step. Computationally, the prob- 
lem addressed in the figure is significantly easier than 
the "open-loop" problem illustrated in Figure 2, as the 
number of active cells by t — 1 is reduced from 50,000 
to only 4,000; the computation time for this simulation 
is also reduced accordingly, from 40 to 4 seconds for the 
time interval shown. The PDF at time t = 1 splits into 
two concentrated regions on the left and right sides of 
the figure. Future measurements might reveal in which 
region the state really is; until such measurements are 
received, the GBEES algorithm is perfectly capable of 
following both. A plain Kalman filter, which assumes a 
central estimate, would necessarily fail to model such a 
splitting. 

5 Analysis and Conclusions 

A novel algorithm is introduced in this paper to ex- 
ploit the remarkable sparsity of the evolving PDF in 
phase space, thereby, for the first time, making high- 
resolution grid-based Bayesian estimation computation- 
ally tractable for nontrivial systems. The method gener- 
alises straightforwardly to any number dimensions, with 
computational cost expected to be a trade-off between 
the curse of dimensionality and the increased sparseness 
of the PDF. In application, the algorithm developed is 
shown to track, with unprecedented fidelity, the com- 
pletely non-Gaussian PDF of the estimate of a Lorenz 
system, both with and without measurement updates. 
The simulation exhibits a competition between infor- 
mation loss due to the random state disturbances and 
stretching of the PDF in the unstable directions of the 
system, and information gain from measurements. 

Grid-based Bayesian estimation algorithms are some- 
times referred to as approximate grid-based methods. 
We point out that the numerical analysis of §3.1 estab- 
lishes that the numerical method used to propagate the 
Kolmogorov equation in the present grid-based estima- 
tion algorithm is second-order accurate in both space and 
time; this compares favorably to the (slower than lin- 
ear) 0(1/ y/N) convergence rate of particle methods ap- 
plied to the Kolmogorov equation (see Bernard, Talay, 

6 Tubaro 1994). 

Finally, Lagrangian (that is, particle-based) simulation 
techniques have been explored for decades in the field 
of fluid mechanics, but for n > 2 remain mostly a re- 
search novelty. On the other hand, grid-based methods 
(often with adaptive grids to focus the computational ef- 
fort where it is needed) have proven immensely success- 
ful in a variety of complex situations in fluid mechanics, 
such as in the characterization of fluid turbulence and in 
the design of commercial airliners, where computational 
methods have largely supplanted repetitive wind-tunnel 
testing. There appears to be no reason why the same 



success of grid-based methods will not also be realized 
in Bayesian estimation approaches, once the community 
working on such problems fully appreciate how the re- 
markable sparsity of the PDF in such problems may be 
exploited. 
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